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Abstract 

In this paper, wc develop a higher order symmetric partitioned Runge-Kutta 
method for a coupled system of differential equations on Lie groups. We 
start with a discussion on partitioned Runge-Kutta methods on Lie groups 
of arbitrary order. As symmetry is not met for higher orders, we general- 
ize the method to a symmetric partitioned Runge-Kutta (SPRK) scheme. 
Furthermore, we derive a set of coefficients for convergence order 4. The 
SPRK integration method can be used, for example, in simulations of quan- 
tum field theories. Finally, we compare the new SPRK scheme numerically 
with the Stormer-Verlet scheme, one of the state-of-the-art schemes used in 
this subject. 
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1. Introduction 

In the simulation of gauge theories in lattice Quantum Chromodynamics 
(QCD), for example, one is interested in calculating expectation values of 
certain operators. That means, very high dimensional (10^ or more) inte- 
grals have to be evaluated. As this can not be done analytically in general. 
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numerical methods are applied to get approximations to these integrals, i. e., 
the expectation values. 

The Hybrid Monte Carlo (HMC) method [3l H] is widely used here. Alike 
in Monte Carlo integration, in HMC the integration is realized by averaging 
over evaluations of the integrand at certain, suitably chosen (importance 
sampling), values. In HMC, these points - or samples - are drawn from a 
combination of so-called Molecular Dynamics (MD) and Metropolis steps. In 
the former, starting from a suitable sample, a candidate for a next sample is 
derived from solving some differential equations. In the latter, it is checked 
whether this candidate suits or not, i. e., if it follows a certain distribution 
or not. 

In this paper, we take a close look at the numerical integration of the 
differential equations in the MD step, arising in QCD problems. Commonly 
the Leapfrog (Stoermer-Verlet) scheme, Omelyan methods p!T| [Ti] or split- 
ting methods with multiple timescales a la Sexton- Weingarten [13] are used. 
We formulate time-reversible higher order integrators that are based on im- 
plicit partitioned Runge-Kutta schemes and show that they allow for larger 
step-sizes than the Leapfrog method. 

The paper is organized as follows: The equations of motion within the 
MD step of HMC, together with the requirements to preserve the Lie group 
structure and time reversibility in numerical integration schemes, are intro- 
duced in section 2. Partitioned Runge-Kutta (PRK) schemes are discussed 
in Section 3, based on Magnus expansion and Munthe-Kaas approach. They 
define numerical methods which preserve the Lie group structure. As this 
class of methods allows only for order lower or equal two, if in addition time 
reversibility has to be met, we generalize PRK methods to symmetric PRK 
(SPRK) methods which allow for higher order and derive a set of coefficients 
for a method of order 4 with 3 stages. The numerical results obtained in 
Section 5 show the efficiency of SPRK methods compared to the Leapfrog 
scheme at lower tolerances. It turns out that the integration measure (area) is 
not preserved with SPRK methods, i.e. they are not symplectic. This means 
that the determinant of the Jacobian has to be included in the Metropolis 
step. A conclusion and outlook to open question and future work is finally 
given in Section 6. 
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2. Equations of Motion 

We will not go into details for the origin of the differential equations to 
be solved in the MD step. For a detailed discussion we refer to [12] and [8]. 
The dynamical system to be solved can be thought of as equations of motion, 
derived from some Hamiltonian operator H{Y, In lattice QCD especially, 
these equations of motion form coupled systems of matrix differential equa- 
tions of the form 

y.^ (la) 

dH (Y 

^. = - ; ' ^ =9uiY), for^. = l,...,n. (lb) 

Thereby, F is a vector of n elements yi, . . . , i/n, each being an element of 
a matrix Lie group G; the vector \l/ comprises n elements ipi, {u = 1, . . . ,n), 
each being an element of the Lie algebra q associated to the Lie group G. 

The coupled system ([T| becomes an initial value problem (IVP) by pre- 
scribing initial values: yu{0) := y^fl ^ G and ipu{0) '■= ipufi G for z/ = 
1, . . . ,n. 

This IVP is usually solved by a numerical integration method on a time- 
grid {to = 0,ti,t2 . . . , W}- Here we let {Yi, "ifi) = ^hiXi-i, represent 
a one-step method that computes an approximation (Yi, \E'i) ^ (Xiti), ^^(t/)) 
at a time-point ti, from values (Yi^i, \E'i_i) of the Lie group and Lie algebra 
elements at the time-point ti-i. The progress in time is given by the step- 
size h = ti — ti-i- The accuracy of the method $/i is measured by the 
deviation ei = ||(y;,\&i) — (Y(ti),'if(ti)\\ with some suitable norm || ■ ||. The 
method is said to be of local order p if e/ = 0{hP~^^) and (y/_i, \I'i_i) = 
iY{U.,),^{U.,)). 

For efficiency, the order p of the method should be preferably high, as 
this allows for large step-sizes h to satisfy prescribed error tolerances. In 
addition, numerical schemes applied to the Lie group problem ([T]) have to be 
equipped with the following properties: 

• The Lie group structure has to be preserved. That means, approxi- 
mations to Y and \E' have to reside in G and q, element by element, 
respectively; 

• The integration scheme has to be symmetric. This is a consequence 
of the detailed balance condition of the Markov process defined by the 
HMC method (see [3]); 
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• Detailed balance requires area preservation, i.e. a symplectic integra- 
tor. If the area is not preserved by the integration, the determinant of 
the Jacobian has to be included in the Metropolis accept-reject step. 

In lattice gauge theories, the state-of-the-art integration methods applied to 
([T]) are the Leapfrog scheme of order 2 as well as Omelyan pT| or splitting 
methods for higher orders |13]. In this paper, we develop an alternative to 
these schemes on the basis of partitioned implicit Runge-Kutta methods. 



3. Runge-Kutta Methods for Lie Group Problems 

Applying a numerical integration method directly to the coupled system 
([T]), it can not be guaranteed that the approximations to (for u = 1, . . . , n) 
are elements of the matrix Lie group G, which is closed under matrix mul- 
tiplication but not under summation. Hence, measures have to be taken to 
preserve the Lie group structure in the numerical approximation. 

3.1. The Magnus expansion 

The theorem of Magnus [9] allows to transform the Lie group differential 
equation (la) to 

ujy = deyo^'li^ilj^), for z/ = 1, . . . ,n, (2) 

with UJy{t) G Q and UJy{^) := 0, i. e., to a differential equation in the cor- 
responding Lie algebra q. The way back from the Lie algebra to the Lie 
group G is given by the mapping Uuit) = exp (uuit)) yufl where exp(-) is the 
matrix exponential. 

The central point in the transformation (|2| is dexp~^, the derivative of 
the inverse of the matrix exponential. This is given by the series 

o^expjj i^y) = (^^) • (3) 

fc>0 

Here, Bk is the k-th Bernoulli number and ad^^ is the adjoint operator, de- 
fined by &d^,{ipu) ■= [uJu,il^u] = ujui^u - i^ui^u and adf^^(V'i.) = [uu,a.d^J^^{iiu)] 
with the convention ad°^('?/'i/) = ipv For a detailed discussion we refer to [5]. 
In total, we record that the problem IVP ([T]) is equivalent to 

fc=0 

Tpy = gy{Y) with Y = {yu),y=i,...,n where y^ = exp(a;^)?/^,o, (4b) 
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with y^{0) := y^fi G G, tp^iO) := tp^fi G and u^{0) := G for z/ = 1, . . . , n. 

This transformed problem can now be solved directly by a Runge-Kutta 
method without destroying the Lie group structure: as the Lie algebra is a 
vector space pLj both analytic solution u^{t),ip^{t) as well as approximations 
attained by a numerical integration scheme are elements of the Lie algebra 
Q. Furthermore, as for any a G the matrix exponential exp(a) is in the 
associated matrix Lie group G, also is in G. 

3.2. The Munthe-Kaas approach 



Clearly, in practical computations the series in (4a), which is the expan- 
sion of (iexp~J given in equation (|3|, can not be evaluated. Instead, one has 
to truncate the series after some q + 1 terms, i. e., sum up for = to A; = q. 

Munthe-Kaas explains in [10] how the truncation index q can be chosen 
properly. According to the observations made therein, a numerical integra- 
tion method of local order p applied to the Lie algebra problem Q, demands 
to take into account at least the first p — 1 addends in the series. That means, 
the truncation index q has to satisfy q > p ~ 2. 



By truncating the series in (4a) at /c = g = p — 2, for a fixed p G N, the 
coupled system ^ is formally replaced by the truncated IVP model 

k=0 

i)^ = g^{Y) with Y = {yu)u=i,...,n where y^ = exp(w^)y^,o, (5b) 



with y^{0) := y^^ G G, V^^(O) := tpifl G and u)^{0) := G for z/ = 1, . . . , n. 

We can now apply a numerical integration scheme of local order p. This 
yields approximations {yiy,i,uju,i,ipiy,i) to the exact solution {y^,Uu,4'u){h) of 
the truncated model ([s]) at time-point ti = + h of local order p. That means, 
the approximations satisfy \\ipu{h) — ipu,i\\ = C(/i^^^), correspondingly for y^, 
and Qu, for z/ = 1, . . . , n. 

The central statement of Munthe-Kaas [lO] is that {yu^iyOJu^y'ipu,!) is also 
an approximation of local order p to the exact solution {y^,co^,4'u)ih) of 
the original problem (|4]). Or, the other way around: a method to compute 
approximations of order p to the exact solution of the Lie algebra problem 
^ consists in applying a numerical integration scheme of order p to the 
truncated model ([s]). 
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In the latter problem, i. e., in the dynamical system (|5|, we skip the 
^-notation and use in the following the condensed formulation 

n = f{n,^), (6a) 

4/ = g{Y) with Y = exp(fi)Fo, (6b) 

where 

Y{t) = (y^(t))^=l,...,n, ^(t) = iMt))^=l,...,n, andfi(t) = Mt))u=l,...,n, 

and 



/(^],vl>) := {UM.)),=,_n with /.(a;.,V^.) = $^^ad^^(^, 
9{y) ■= {9v{y))^=i,„n and exp(fi)Fo := (exp(u;^)?/^,o)i=i ... „ 



with initial values F(0) := Fq := {yvfl)v=i,...,n, ^(0) := := (^i.,o)v=i,...,n 
and fio := (0)i.=i,...,n- 

3.3. Partitioned Runge-Kutta methods 

To solve the coupled dynamical system (|6]), we apply a partitioned Runge- 
Kutta (PRK) method [6] with s stages and coefficients hi, hi, aij, cxij for i, j = 
1, . . . , s. Starting from t = we first compute approximations Vli and at 
the time-point t = hhy 

s s 

ni = no + hY^ b^Ki, = ^'o + h^hiLi, (7a) 

i=l 1=1 

with increments Ki and Li for i = 1, . . . , s defined by 
where Oj, ^i and 1^ are internal stages given by 

s s 

Vti = VtQ + h^^ aijKj, = \[/o + /i ^ (^ijLj, % = exp (Hj) Fq- (7c) 
i=i i=i 
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Finally, we get an approximation {yu,i,'4'u,i) of local order p to the exact 
solution {uuit), Tpyit)), one time-step h ahead, i. e., at t = + /i by 

yu,i = (l^i), with Yi := exp (fii) Fq and = (^i)^ . (7d) 

The coefficients 6j, foj, aij, ctij for i, j = 1, . . . , s steer the behavior of the 
method and have to be chosen properly. Conditions, the coefficients have 
to satisfy to attain an approximation of local order p, are found by series 
expansions of the approximations produced by the method ([T]) and the exact 
solution of the truncated problem ([s]) in powers of h, followed by a comparison 
of the series' coefficients. For an accuracy of order 2, for example, the order 
conditions for p = 1 and p = 2 as stated in Tab. [TJhave to be fulfilled. 



p 




m 


1 




= 1 


2 


E^b^ai = 1/2 


= 1/2 



Table 1: Order conditions for the PRK method ([t]), see [S] or [B]. 



4. Symmetric PRK methods for Lie group problems 

Symmetry, frequently termed time-reversibility, is closely related to the 
adjoint of a method. The adjoint method $^ of a method $/j is the inverse 
of the original method with reversed time-step —h, i. e. 

'^>*. := ^Zl (8) 

A numerical integration scheme $/i is symmetric if it equal to its adjoint 
i. e., $/i = The specifications for the adjoint method are found by 
exchanging the roles of the initial value and the approximation, reversing 
the time-step and dissolving this system again for the approximation. For a 
detailed discussion of the procedure, we refer to [5]. 
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4-1- Symmetric Lie group PRK method 

When deriving the adjoint method of the PRK method ([T]), however, we 
have to pay attention to the Lie group structure. Therefore, we state some 
more details. 

As described, we exchange (/i,lo, ^0,^1,^1) in ^ by (-/i,Fi,^i,Fo,^o) 
along with replacing Vti with —Vti and then resolve for (Yi, \E'i). 

The arising instructions can be stated as a PRK method, similar to the 
original method ([T]). The coefficients 6*, 6*, a*-, a*- for i,j = l,...,s are 
connected to the coefficients of the original method via 

bi '■= hs+i-i, and a*- := bg+i-j — Q;s+i-i,s+i-j, 

and for 6j and aij in equal measure. There are, however, two important 
differences. The increments Ki in the adjoint method are given as 

s 

Ki = /(a, ^^), with Ui := h J2i-(^s+i-i,s+i-j)Kj. (9a) 

^* is defined in the obvious way. Furthermore, the internal stages for the 
links are defined by 

Y,* = exp{n,)-exp{n,)Yo. (9b) 

The numerical scheme ([T]) is symmetric if it coincides with its adjoint 
method. That means, the approximations {Yi, "$1,0,1) produced and, there- 
fore, the increments [Ki, Li) within $/i and have to agree. 

Immediately, the conditions 

bi = bs+i-i, (10a) 
bi = bs+i-i, (10b) 
aij = bs+i-j - as+i-i,s+i~j (10c) 

for dAl i,j = 1, . . . ,n become clear. Considering the increments Li, we see 
that ^i = ^i needs to hold. However, we can only derive a condition for 



the coefficients if fli and fii in (9b) commutate. The condition imposed then 
reads 

aij = bs+i-j - as+i-t,s+i_j. (lOd) 
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fli and fli commutate if is some multiple of —ag+i^i^s+i-j, i- e., if 



bj = di ■ aij with some di G 



(lOe) 



Finally, we consider the increments Ki. Comparing (9a) and (7b) yields the 
condition 



a,; 



'0:s+l-i,s+l-j- 



(lOf) 



The conditions (lOa-d) are the usual symmetry conditions for PRK methods 



(see [5]). The conditions (lOe f) arise from the Lie group nature of the 
problem. 

The symmetry conditions (lOd) and (lOf ) imply 6j = for alH = 1, . . . , n. 



Clearly, this is a contradiction to the basic order- 1 condition fej = 1. 

However, this contradiction does not appear if the local order p of the 
method is 2 at most. In this case, (5a) simplifies to Ui, = ip^,. Hence, the 



function / only depends on i. e., / = f{^)- Consequently, the symmetry 



condition (lOf) vanishes, such that the conflict between order and symmetry 
condition disappears. 

This statement is summarized in the following 

Lemma 4.1. The PRK scheme ^ can only be symmetric if p < 2. 

Indeed, there is a symmetric partitioned Runge-Kutta scheme for Lie 
group problems. The Stormer- Verlet scheme, also known as Leapfrog method 
in various applications, can be interpreted as partitioned Runge-Kutta scheme 
of type ([T]) with coefficients given by the extended Butcher tableaus (see [5]) 



and 



a 


A 




b 












1 


0.5 


0.5 




0.5 


0.5 



a 


A 




b 



0.5 


0.5 





0.5 


0.5 







0.5 


0.5, 



where A = (ajj)i,j=i,...,s, b = {bi,...,bs) and a = (ai, . . . , a^)^ with 
tti = aij, and in the same way for A, b, a. 

It is easily checked that these coefficients define a method of order two 
applied to the original system and satisfy the conditions (10) with di = 
and = 0.5 in equation (lOe). Therefore, the Stormer- Verlet scheme is a 
symmetric integrator. 
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4-2. Higher Order Symmetric Partitioned RK Methods 

To allow for higher order symmetric methods, we need more flexibility in 
the stage vector Yi. We get this by introducing additional coefficients for 



j = 1, . . . , s and replacing the links' internal stages Yi in (7c) with 
Yi = exp (Xi) Yo = {exp{x^^i) ■ Z/i.,o)^=i,...,„ , 

s 

where X, = 



(11) 



By this, the aforementioned contradiction in the conditions vanishes be- 



cause (lOf) and (lOe) are replaced in a first step by 



— bs+l-j — 7s+l-i,s+l-i 

and bj = di--fij. 



(12a) 
(12b) 



However, deriving the order conditions for this adapted PRK scheme. 



it turns out that the condition (12b) for symmetry leads to a conflict with 



conditions for convergence orders p > 2. The problem is caused by the term 
exp(ni) that enters the definition of the internal stage variable Y^ of the 
adjoint method. The coefficients bj have to be the multiple of some other 



method's coefficients to guarantee the commutativity of Qi with Qi in (9b) 



or with Xi = h^A—'ys+i-i,s+i-j)Kj if the variant (11) is used. 



As a solution, this necessity - and therefore the condition (lOe) and (12b), 



respectively - can be avoided by a further reformulation of the stage vector 



Yi. Finally, we replace the link's internal stages in (7c) with 



Yi = exp (Xi) exp (|fii) Yq = {exp{xu,i) ■ exp(cu,,,i) ■ 2/z.,o)^=i^,.._„ , 

s 

where Xi = {xu,i)^=i^_^n = ^^^Ihi^i- 

and arrive at the symmetric partitioned Runge-Kutta scheme (SPRK) 



ni = hY,biKi, ^^ = ^o + hY,biLi, Fi = exp(l]i)ro (13a) 



j=i 



i=l 
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with increments Ki and Lj for i = 1, . . . , s defined by 



and internal stages defined as 



(13b) 



Yi = exp (Xi) exp (^l^i) Yq, = /i ^ Tij^j, 



(13c) 
(13d) 



with coefficients bi, bi, aij, Sij, jij for i, j = 1, . . . , s. 

The SPRK method (13) is not a "standard" partitioned RK method. 



Hence, we have to compute both symmetry and order conditions for this 
method. 

4-2.1. Symmetry conditions 

We first determine the conditions for symmetry that the method's co- 
efficients have to fulfilL For this, we determine the adjoint method of the 



method (13) and carry out the same steps we described in Sec. 4.1 for the 



method at hand. 

Again, the adjoint to the SPRK method can be formulated as a parti- 
tioned RK method with coefficients bi,bi,aij,aij (like in the adjoint for the 
PRK method ([t])) plus coefficients 7*- := —'ys+i-i,s+i~j for i, j = 1, . . . , s. 

We recall that the root for the problems we recorded was the computation 
of the PRK adjoint method's increments Ki and the internal stages Y* as 



stated in (9a) and (9b), respectively. Within the adjoint method of the SPRK 



scheme, the Kj's are computed in the same way (9a) but the internal stages 
Y* now amount to 



F,* = exp {Xt)-exp{-ln,)-Y^ 
= exp (X,*) ■exp(ifii) -Fo, 



(14) 



with X* defined in the obvious way (as hj^jltj^j). 

Finally, the symmetry conditions, which the method's coefficients have 



11 



to satisfy, are 



bi = bs+i^i and bi = b^+i-i, 

(15) 



Oiij — bs+l-j — as+l-i,s+l-jy 
lij = —ls+l-i,s+l-j 

for all i,j = l,...,s. 

Note that the absence of any condition like (|10e ) or (12b) is due to sim- 



ilar form of the specifications of Yi and Y* in (13d) and (14), respectively. 
Because of that, ^Qi does not have to commutate, neither with Xi nor with 

Furthermore, the symmetry conditions do not conflict with the order 
conditions that we derive now. 

4-2.2. Order conditions 

As explained by Munthe-Kaas, the SPRK method ([13]) is of local order p 
with respect to the Lie group differential equation (jl]) (see Sec. [3]) if 

\\Y{to + h)-Yi\\=0{hP+^) and ||^(to + Z^) - ^i|| = C(/i''+^), (16) 

where Yi and are the approximations to exact solution Y{tQ + h), \I'(to + /i) 
of the suitably truncated problem ^ after one step. 

As within the SPRK method, the approximation to the links Yi arises 
from evaluating the matrix exponential (we assume here that we can evaluate 
this exactly), which is Lipschitz on every closed interval, it suffices to demand 

\\Q{to + h) - QiW = 0{hP+^) and ||^(to + /i) - ^i|| = C(/i^+^). 



As mentioned in Sec. 3.3[ expanding (f2i, \E^i) and (fi(to + h), \l/(to + h)) to 



power series in h leads to order conditions for the method's coefficients by 
matching the series' coefficients of the exact solution and its corresponding 
approximation. In constructing RK schemes, the use of some B-series' ap- 
proach [Tlin] is common. However, to our knowledge, there is not yet a theory 
available that fits to the situation reported in this paper. Therefore, for the 
time being, we apply standard Taylor series expansion. 

We are concerned with a non-abelian, i. e., non-commutative, matrix 
Lie algebra q. Therefore, the Taylor expansion and the nature of the order 
condition deviate from the standard, i. e., abelian case. We recognize several 
issues: 
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The right hand side (6a), i. e., the function /, depends on the maximal 
order p we aim to attain. Each increase in order adds an additional 
fc-th order commutator ad^^{ilJu)- As a consequence, additional (w. r. t. 
the abelian case) order conditions appear. 



• The definition of the internal stages Yi in (13d), where the approxima- 
tion Qi is included, is a non-standard RK formulation. This and the 
non-commutativity of Xi and fli also leads to further order conditions. 

As the actual computation of the Taylor series's expansion in this case 
is tedious, we skip details. Instead, we just state the order conditions up to 
order p = 3 in Tab. [2j Note that, due the maximum order p = 3, the function 
/ had to be chosen as: 



p 




m 


1 




T.A = i 


= 1 




2 








1/2 


3 


T- ■ 


bi {aiOi - aijOij) = 1/6 




1/3 






«i.(7, + EfcW2) = l/6 


Y.i,fii{lij + bj/2)aj 


= 1/6 



with 7i := Y.. 'jij and := Y.j 



Table 2: Order conditions for the symmetric partitioned Rmige-Kutta (SPRK) method 
which is stated in (13). 



The following set of coefficients solves both the order conditions given in 



Tab. 2 and the symmetry conditions specified in (15) and, hence, defines a 



13 



SPRK method of order p = 3 with s = 3 stages: 

«21 = 7ll = -V3/6, 023 = 733 = VS/Q, &2 = 1, 

Sii = (3 + v^)/6, a2i = (3 + v^)/12, 61=63 = 1/2, 
a23 = (3 - v^)/12, S31 = 1/2, 



"33 



-V3/6. 



Here, coefficients that are not mentioned exphcitly are 0. It is known in the 
hterature [5J that symmetric Runge-Kutta schemes of order higher than 2 
as well as partitioned ones for problems of the type ^ can not be explicit. 



Indeed, the coefficients (17) define an implicit SPRK method. It follows that 
the equations ([l3|3-d) for the increments Ki and Lj of the SPRK scheme are 
implicit. We solve them through a fixed-point iteration. 

4-3. Global error 



Recall that the SPRK method (13) is a numerical scheme that is con- 



structed such that it produces approximations of local order p with respect 
to the exact solution of the truncated IVP But, due to the argumen- 



tation of Munthe-Kaas (see Sec. 3.2), the results are also approximation of 



local order p w. r. t. the original problem, i. e., the equations of motion ([T|. 
That means, (16) holds also for the exact solution (Y(tQ + h), ^(to + h)) of 

&■ 

We now investigate the global error, i. e., the deviation between the exact 
solution (Yito + N ■ h), "^(to + N ■ h)) and the approximation (y^r, \l/7v) after 
steps of the numerical scheme with step size h. 



It is readily checked that the SPRK method ( 13 ) is consistent in the sense 
of Theorem 8.1 in [6l II. 8] and, hence, has an asymptotic expansion of the 
form 

\\Y{tN)-Yr,\\ = 

evAtNW + ey,p+i(t,v)/i^+' + ■ ■ ■ + ey,fc(t^)/i^ + EY,h{Y,tN)h^^\ (18) 

Thereby, the eY,v{tN) are solutions of some differential equations and Ey^h is 
bounded on [to, tend] (accordingly for \E'). Furthermore, it holds = to+N-h. 
As, in addition the SPRK method is symmetric. Theorem 8.10 in [6, II. 8] 



applies, which says, that the expansion ( 18 ) only contains even powers of h: 



\\Y{tN) - YnW = eY,2,{tN)h'' + ey,2,+2(tjv)/i'^+' + • • 
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In view of the third order SPRK method with 3 stages defined by the set 
of coefficients (17), we therefore have a global error 



where \Ef(t)) is the exact solution of the equations of motio n ^ and 

{Yj\f, '^n) are the approximations produced by the SPRK method (|13[). 



5. Numerical Results 

We test the SPRK method of order 4 numerically. For this purpose, we 
simulate an SU(2, C) gauge field, used in lattice Quantum Chromodynamics, 
by means of the Hybrid Monte Carlo method. The equations of motion ([T]) 
are generated by an Hamiltonian 

HiY,^) = E^^^i^) + SGiY) 

with kinetic energy Ey^^^ and so-called Wilson action Sg{Y). This is in 
detail described in |2l paragraph 7.2.3]. Thereby, Y = {yu)u=i,...,n is a field 
of elements of the special unitary Lie group SU(2, C), and \E' = {4'u)i'=i,...,n a 
field of elements of its associated special unitary Lie algebra su(2,C). 

In lattice gauge theories, the elements of Y are called links. Each link y^, 
has an associated fictitious momentum which is a traceless and hermitian 
2-by-2 matrix. Thus, the momentum is connected to traceless and ant- 
hermitian element ipu of the Lie algebra 5u(2,C) via a multiplication with 
the complex i: = {ip^ = ipu)u=i,...,n ■ 



For completeness, note that the function QuiY) from equation (lb) is not 
just evaluated at the lattice point v itself but also on some adjacent lattice 
points called staples. However, this fact is not important for the derivation 
of the SPRK method. 



The simulation is performed on a 2-dimensional lattice of dimension 8x8. 
As energy is preserved analytically along the trajectory, (|Aif|), the mean of 
the absolute difference between the numerical approximation to the Hamilto- 
nian at the end of each trajectory of unit length 1 and the initial Hamiltonian, 
is an easy to get measure for the numerical approximation error. Figure [T] 
reveals a global error 0{h'^) of the SPRK method compared to an error of 
order 0{h?) for the Leapfrog scheme. This fact is reflected in Figure [2| show- 
ing that the SPRK method is more efficient than Leapfrog with respect to 
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CPU time at smaller energy violations {\AH\). For lower error tolerances, 
however, the numerical effort for fixed-point iteration in the SPRK method 
becomes visible, and the efficiency is drastically reduced. The area preserva- 
tion mentioned in section [2] is as expected not met: numerically we have an 
error of order 4 for the determinant of the Jacobian d{Yi^, fyf)/d{Yo, \l/o)- 
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Figure 1: Numerical approximation error of (|A77|) for Leapfrog (+) and SPRK (x) 
scheme for different step sizes. The mean of the energy change |AiJ| along a trajectory 
with length 1 is computed from a simulation that is comprised of 5000 trajectories on a 
2-dimensional 8x8 lattice. 



10" 



S 10- 



I 10"' 
O 



10"' 



'X, 



+ Leapfrog 
X SPRK4 



10" 



10" 



X + 



XX 



x 



10"= 
<|AH|> 



10" 



10" 



Figure 2: CPU time versus accuracy for Leapfrog (+) and SPRK (x). These values are 
measured in the aforementioned simulation on a 2-dimensional lattice of the size 8x8. 
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6. Conclusion and Outlook 



In this paper, we have developed symmetric partitioned Runge-Kutta 
schemes for matrix differential equations of type ([T]), which preserve both 
Lie-group structure and time-reversibility, and allow for higher order at 
the same time. Especially for smaller energy violations (higher error toler- 
ances), SPRK schemes have turned out to be more efficient than the Leapfrog 
scheme. However, for larger energy violations (lower error tolerances), the 
fixed-point iteration turned out to be the computational bottleneck. This 
drawback has to be attacked in future works. One idea is to replace the 
somehow artificial term exp{Qi/2) introduced to allow for symmetry, which 
causes a strong coupling of all components and thus is responsible for a high 
degree of non-linearity in system (13). As symplecticity is another desirable 
property, which simplifies the acceptance step within HMC, the derivation 
of symplectic SPRK schemes is a next natural step. 
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